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Abstract. Dark matter halos are the building blocks of the universe as they host galaxies 
and clusters. The knowledge of the clustering properties of halos is therefore essential for 
the understanding of the galaxy statistical properties. We derive an effective halo Boltzmann 
equation which can be used to describe the halo clustering statistics. In particular, we show 
how the halo Boltzmann equation encodes a statistically biased gravitational force which 
generates a bias in the peculiar velocities of virialized halos with respect to the underlying 
dark matter, as recently observed in N-body simulations. 
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1 Introduction 

Over the last decade, we have accumulated a good deal of observational evidence that the 
large scale structure (LSS) originated from seed fluctuations in the very early universe, which 
grew via gravitational instability at later times. In this picture, galaxies form within Dark 
Matter (DM) halos [1], and the goal of the theory of bias [2] is to relate observable properties 
of tracers such as the density contrast of galaxies to the underlying matter distribution and, 
ultimately, to the initial conditions. 

DM halos are, to a large extent, the building blocks of what we observe. Knowing 
how DM halos form and evolve under the action of gravity is crucial to understand the 
environment that harbors galaxies and clusters. Understanding the clustering properties 
of halos represents, therefore, a key ingredient towards an accurate description of galaxy 
clustering statistics. 

The aim of this paper is to write down an Effective Boltzmann Equation (EBE) for 
the single phase space density of halo centers, which can be used to derive ensemble average 
halo clustering statistics. Therefore, this must be an effective description, in which the halo 
distribution should be thought of as being the mean-field distribution associated with a given 
realization of the DM distribution. This is the correct interpretation of the halo overabun¬ 
dance defined through a perturbative bias expansion so long as the surveyed volume remains 
finite, as is the case of real or simulated data [3, 4], 

Since the halo Boltzmann equation will describe the dynamics of the halo mean-field, 
one should expect to recover exactly the Boltzmann equation for DM only if the halo centers 
are statistically unbiased relative to the coarse-grained DM density field. This is assumed 
to be true in the coupled-fluids approximation for the coevolution of DM and halos, which 
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is largely used to compute the time evolution of bias by following the evolution over cosmic 
time and in Eulerian space of the halo progenitors - the so-called proto-halos - until their 
virialization [5-8]. In this framework, the momentum equations for halos and DM do not 
differ at all. However, this approximation is not valid in general, as one should expect that 
the force acting on large scale structure tracers is indeed biased relative to the force acting 
on a random held point. In other words, the sampling bias induced by the tracers propagates 
into their dynamical evolution. 

Therefore, the EBE for halos will be different from that of DM. This point should not be 
underestimated: the DM Boltzmann equation is the starting point of any analytical approach 
to the description of the LSS (see e.g. [9]). Here we advocate that, since the ultimate goal 
is the understanding of the statistical properties of what we observe (that is galaxies), one 
ought to adopt the EBE discussed below as the starting point for halo centers. 

This is not the only reason for pursuing this project, though. LSS probes such as 
peculiar velocities or redshift space distortions can be used to test the law of gravity through 
a measurement of peculiar velocities, i.e. deviations from pure Hubble flow [10, 11]. However, 
peculiar velocities are effectively measured at the position of halo centers. Therefore, even 
if the latter locally flow with the DM, velocity statistics may be biased if the halo peculiar 
velocities do not provide a fair sample of the matter flows (for recent discussions see Refs. 
[12-17]). Thus far, measurements of the halo velocity power spectrum appear consistent with 
little or no statistical velocity bias [13, 18]. However, [12, 19] reached somewhat different 
conclusions. Furthermore, the effect must vanish in the limit k —> 0 so that, at low redshift, 
scale-dependent nonlinearities complicate the interpretation of the results. 

Peak theory predicts the existence of such a statistical effect which, in linear theory, 
manifests itself as a reduced halo velocity dispersion [20] together with a ^-dependent velocity 
bias in the usual Kaiser formula [21] and in the time evolution of the linear peak bias [3] . 
The authors of Ref. [19] have recently detected through N-body simulations a large scale halo 
velocity bias b v (k ) 


b v (k) = (l-R 2 v k 2 ), (1.1) 

which appears to remain constant throughout time until virialization. Here, R v is the typical 
scale of the halo velocity bias 

Rl = a] = J -0^ k 2 ip(k)W\kR), (1.2) 

P(k ) is the DM power spectrum and W ( x ) is a spherically symmetric, time-independent kernel 
with filtering scale R equal to the Lagrangian halo radius. This result confirms that, when 
the purpose is to compute statistical quantities such as power spectra, etc, the coupled-fluid 
approximation requires a modification of the momentum conservation equation for halos [19] 
as the force felt by the halos is biased relative to the force acting on a random DM particle. 

The EBE provides a way of describing the evolution equation for biased tracers. As 
first suggested in Ref. [22], it encodes the statistical velocity bias which originates from the 
fluctuations around the ensemble average of the gravitational force acting on the halo and it 
predicts the right value (1.1) observed at the linear level. Therefore, the EBE explains at more 
fundamental level why and how one needs to modify the momentum fluid equation for the 
coupled-fluids of halos and DM, as suggested by Ref. [19]. In addition, it allows us to compute 
the halo velocity bias at higher orders and to derive additional stochastic components that 
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describe in a mean-field way (i.e. at the statistical level) how biasing influence the dynamics 
of halo centers. 

The paper is organised as follows: in Sec. §2, we motivate our approach and introduce 
the various averaging we perform throughout the calculation. In Sec. §4 we perform a detailed 
computation of the left-hand side of the EBE adopting a general method based on the path 
integral approach for computing conditional averages. Section §5 is devoted to a similar 
discussion regarding the right-hand side of the EBE. We then conclude with some final remarks 
in Sec. §6. 


2 From the Klimontovich-Dupree equation to the Boltzmann equation 


Let us first explain in more details what we mean by halo “mean-held”. In a local bias 
approximation (e.g. [23]), one can express density fluctuations in the halo mean-held as the 
series expansion 

5 h (r,t) = bi5 R {r,t) + ^b 2 (S R (f,t) - al(t)) + ... (2.1) 

where the filtering scale R is proportional to the halo mass [21], the bias parameters b±, b 2 
etc. are ensemble average whereas 8 R (r,t) is one particular realisation of density fluctuations 
in the large scale structure at cosmic time t. The filtering rehects the fact that we are not 
interested in the internal structure of the halos, only in their center-of-mass position and 
kinematics. From the knowledge of this halo mean-held, we can compute halo clustering 
statistics, such as the halo 2-point correlation, as 

6i(®i - x 2 ,t) = bl(8 R (x 1 ,t)5 R (x 2 , t)) + ... (2.2) 

where (...) is an average over realisations of the large scale structure. 

Our goal is to write down an EBE for the single particle phase space density /h of halo 
centers, such that its zeroth moment yields 5h(r,t). This EBE must be different from the 
Boltzmann equation of DM: the latter is recovered only if the halos are statistically unbiased 
relative to the coarse-grained DM density held (which they are not). Our starting point is 
the Klimontovich density [24] 


f K ( r , v , t ) = [ r - fi ( t )\ <5 d [ v - Vi ( t )] , (2.3) 

i 


which is the single particle phase space density for DM. We have assumed that all the ele¬ 
mentary DM particles have equal mass m, and we normalize our units so that m = 1. 

The corresponding equation of motion for such a phase space density is the Klimontovich- 
Dupree equation 


<9/k 

dt 


+ v- 


<9/k 

dr 


V$ K • 


d/k 

dv 


= 0 , 


(2.4) 


where 

V4> k (d t) = G n f d 3 r'd 3 v'f K (r',v',t)^— J '\ 
J r — r'\ 


(2.5) 


is the corresponding gravitational force. At this stage we have not learned much, as the 
Klimontovich density encodes the trajectories of all the particles in the system, which is 
much more information than we can handle and want to keep track of (in particular, we are 
not interested in the internal properties of virialized structures). Moreover, eq. (2.4) describes 
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one particular realization of a universe, which may not resemble our own Universe. Therefore, 
we go one step further and average over realizations of universes that have the same large-scale 
structure (i.e. similar coarse-grained properties). In this way we can construct macrostates 
by averaging over a statistical ensemble of microstates with similar phase space density in 
(small) volumes containing a sufficient amount of particles. 

Let us explain in more detail how this ensemble average is performed with the help of fig. 1. 


Microstates 


(a) f K {r,v,t) 





( 

smearing 

tb 


m * 

*** 

(b) Pf(p, v, t) 




Hp, v, t) 


Macrostate 


(c) f{f,v,t ) 



Figure 1. A pictorial exemplification of the ensemble average described in the text. We consider all 
the configurations ( microstates ) for the distribution of DM particles that yield to indistinguishable 
large scale structure configurations (a given macrostate), (a) The distribution of DM particles is 
described in the phase space by the Klimontovich density /K(r, i7, t). Four different realizations are 
shown, (b) After a suitable smearing (the convolution with a window function) we can describe 
these microstates in terms of a smooth density p. In any point r, p can vary in a range of values 
among the microstates: this stocasticity is encoded in the probability distribution P?(p, v, t). (c) 
The macrostate is described by the phase space density f(r,v,t), which follows from the ensemble of 
microstates through the ensemble average (/k(d v, f)) m icro) or equivalently by computing the expected 
value of p with respect to the probability density P?(p,v,t). 

The Klimontovich density describes a set of points in the phase space. Let us focus 
on a small box in the configuration space. In the part (a) of fig. 1 we show four different 
realizations of this box, which yield the same LSS. Their ensemble average, which we denote 
by (• • •)micro) gives the mean value f(r. v, t ): 

(fK(r,v,t)\ = (= f{r,v,t) , (2.6) 

\ / micro \ ^' / micro 

i 

An alternative way of understanding this average is to smear out the discrete distribution 
of points into a continuous density with an associated probability Pp(p, v, t ), where r is any 
position inside the box, as represented in the part (b) of fig. 1. The range of values assumed 
by p over the ensemble defines the probability density Pp(p, v, t ) for having a velocity v and 
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a density p in a small volume centered at position r. provided that the LSS is the one fixed 
in (c) in fig. 1. 1 In this case, f(r,v,t) is given by the mean value of p with respect to the 
probability distribution P?(p,v,t): 

/*oo /*oo 

f(f,v,t)= d ppPp(p,v,t)= dppP?(p\v,t) P?(v,t) = (p\v) P?(v,t). (2.7) 

Jo Jo 

We will drop the subscript r here and henceforth, as it is clear that P(p,v,t ) is defined at a 
given spatial location. 

The distribution function f(r,v,t ) that characterizes the macrostates defines, upon 
taking moments of the velocity v, the overdensity held, the bulk velocity, etc. of a given 
macrostate. For instance, the matter overdensity held 5(r,t ) is the zeroth moment 


p(l + 6(r,t)) 


1 


d 3 vf(r,v,t). 


( 2 . 8 ) 


This distribution function satishes the equation 


df _ df df -4 -4 

-^7 + v - VT- -4 = Vtf-F, 
at or dv 


where 


and 


V<b = 


G n J d 3 r'd 3 v , f(f',v',t) 


(r — r') 




(2.9) 


( 2 . 10 ) 


F(r, v, t ) 


Cov 


V^K{r,t),f K (r,v,t) 


micro 



= G n / 

d 3 r d 3 w 3 f 2c {r,v,r ,v ,t) , 

( 2 . 11 ) 


J 

r - r'\ 


where 

f 2 c(r,v,r',v',t ) = 

f 2 {r,v,f',v',t) - f(r,v,t)f(r',v',t ) 

( 2 . 12 ) 


is the irreducible two-particle correlation function. 

Equation (2.9) is our starting point towards a derivation of the EBE for halos, as it 
enables us to apply constraints involving <5 etc. to define biased tracers of the large scale 
structure. 


3 Halos as extrema of the smoothed density field 

In order to write down an EBE for the halo phase space density fh(r,v,t), we follow 
the evolution over cosmic time, and in Eulerian space, of the so-called proto-halos until their 
virialization. The proto-halos are the progenitors of isolated DM halos, i.e. DM halos that are 
not contained in any larger halo. While their shape and topology change as a function of time 
(smaller substructures gradually merge to form the final halo), their centre of mass moves 
along a well-defined trajectory determined by the surrounding mass density held. Hence, 
unlike virialized halos that experience merging, by construction proto-halos always preserve 

x Note that this smearing should not be confused with the filtering procedure introduced in section 4. The 
latter defines a UV cutoff in momentum space denoted by A, while the former represents an average over 
microstates in small, real space volumes of size £, with the requirement that £ 1/A. 
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their identity. Their total number is, therefore, conserved over time, such that we can write 
a continuity and Euler equation for their number density and velocity, respectively. 

We will consider a specific model for the halo biasing in order to derive explicit expres¬ 
sions for the amplitude and scale-dependence of the effective corrections to the halo EBE. 
Namely, we will hereafter designate by halos (or peaks) those clustered objects that are lo¬ 
cated at the points where V5/j(r) = 0 and where the smoothed density contrast satisfies 
5n(f) = u(Tq (we will occasionally write <5 p k = va o). The subscript R represents quantities 
that have been smoothed on the (comoving) size R of a halo. For instance, 

-j A y 



S(f',t)w(^\r-f'\^ , 


(3.1) 


where W(x) a spherically symmetric smoothing kernel. Moreover, v is the peak height and 
(Jo = do (t) is the variance of the DM field at time t in a given region of size R. In general, 
the spectral moments of the smoothed fields are given by 

*? = I ^^P{k)W\kR), (3.2) 

where W(x) is the Fourier transform of W{x ), P(k) = (\S(k, t)| 2 ) is the matter power 
spectrum, the ensemble average is over the macrostates and the time dependence is implicit. 
We will occasionally refer to those conditions as the peak constraint, although we enforce the 
extremum condition solely. However, notice that, in the high peak limit v 3> 1, nearly all 
extrema are maxima of 5r. 

Last but not least importantly, the filtering of the density field reflects the fact that halos 
are extended objects. In our model, we will take this into account upon assuming that halo 
centers feel a force which is a smoothed version of the force acting on dark matter particles. 


4 The halo Boltzmann equation: the left-hand side 


Our goal is to specialize the Boltzmann equation (2.9) from the DM phase space density 
to the peak phase space density so as to describe the evolution of halos. For this purpose, 
,/ p k(r, v, t) will designate the peak phase space density - i.e. the number density of peaks of 
velocity v at the position r and time t - such that its zeroth moment corresponds to the halo 
mean-field dh(T, t). The density / p k can be written in terms of conditioned probabilities, as 
we will show in eq. (5.2) of section 5. The peak constraint in r also affects the computation of 
the force on the right-hand side of Boltzmann equation, introducing conditioned probabilities. 
As we will discuss in section 5, this side of the Boltzmann equation eventually vanishes. 

Let us focus now on the left-hand side of (2.9). There, the following operator appears: 

V$(r,i) • dvf(r,v,t ), (4.1) 

which is the product of two operators evaluated at the same point, that is, a composite 
operator. It is well-known that, owing to the local product of fields making up the composite 
operator, new ultraviolet divergences 2 appear and depend on the ultraviolet cut-off A ~ 

2 In the picture in which A represents the smoothing scale at which the peak is defined, A never truly 
diverges in CDM cosmologies since it is a physical scale. However, in the formalism we introduce here we wish 
to appeal to the reader’s intuition of the QFT renormalisation scheme. Within this context, our terminology 
is similar to [25]. 


6 





0(1/R). The appearance of these divergences requires the introduction of new counterterms. 
The physical reason why the term (4.1) changes in the equation for halos is that the force 
felt by peaks will be statistically biased with respect to the DM case. The technical way by 
which we can reconstruct the form of V<3? • d$f for the halo Boltzmann equation is through 
renormalization. 

Renormalized operators can be defined and generically expressed as linear combinations 
of all the bare operators of equal or lower canonical dimensionality (see, for instance, [26]) 

[O i (r,t)] = ^2z ij (A)O j (f,t). (4.2) 

3 

Therefore, any composite operator can be expressed as a sum of operators allowed by the 
symmetries of the problem at hand. The renormalized operator [V4> • d$f] being a scalar 
field, it can mix with all the operators allowed by the extended Galilean symmetry [27] and 
must remain a function of the velocity v and preserve the derivative d Vi . We can classify the 
operators Oi that mix with V4> • d$f according to their number of derivatives in v: 


VT • dtf 


V4* • dtf + Zqi V<5 ■ fMf + ... + ZVT ■$/ + .... 


leading order in d$ 


higher order in djj 


(4.3) 


Once we have accomplished this classification, we select the leading operators according 
to the number of derivatives (we are considering large scales) and fields contained in them 
(adding additional fields brings to a result including more power spectra). The fields we are 
dealing with are 4> (which can only appear with derivatives, since the gravitational potential is 
not observable and does not directly affect the dynamics) and the density contrast <5 = V 2 < h/a, 
where a = 3H 2 Q m /2. We can then organise the Oj according to the number of extra couples 
of derivatives (in order to build a vector, we must add two derivatives to V<f> and contract 
two indices) and extra fields appearing in them. The result is shown in fig. 2, where moving 
to the right means adding two derivatives, and moving downwards means adding a field. The 
appearance of derivative operators makes clear that the statistical bias manifests itself only 
through /c-dependent corrections and, therefore, vanishes in the limit of small momenta. 


add 2 der. 

O 0 = V4> • dfff -* Oi = VS ■ d 9 f 

add a field 

0 2 = 5VT • d 9 f _ > 

o3 = (V$-V)V$ • dtf 


Figure 2. Operators allowed by rotational invariance that can mix with V4>-d^/, classified according 
to their number of fields and derivatives. 

The renormalization procedure requires a suitable prescription. For our purpose, it is 
convenient to adopt the following renormalization condition: the correlators of the renormal- 
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ized quantities should reproduce the correlators at the peak. 


( Pi(r,t)\ ■ [0 2 (f,t)] ■ ■ • [O n (r,t)} ) = (Oi{r,t) ■ 0 2 (r,t) ■ ■ ■ O n (r, t)) . (4.4) 

Here and henceforth, the notation (• • • ) p k will designate conditional averages at the peak po¬ 
sition. In our problem, it is precisely the renormalization of the operator V<L(r, t) ■ d$f(r, v, t ) 
that leads to differences in the clustering properties of halos and DM, as we will now see. 
Notice that, since we want to follow the evolution over cosmic time of the so-called proto-halos 
until their virialization, the renormalization condition has to be imposed at all times. 

4.1 Renormalization at the linear order 

In this subsection, we study the renormalization of the composite operator V4>-d^/(r, v. t) 
at the linear order in perturbation theory. 

We begin with the computation of the average of V4>(r, t) ■ d^f(r,v,t) at the peak 
position: 

(v$(r, t) ■ d$f(r, v, t)) = d$f pk {r, v, t) ■ (v$(r, . (4.5) 

This amounts to computing the correlators using the appropriate conditional probabilities. 

Let us first operate at the linear level in perturbation theory such that, at early times, the 
DM density 5 and the gravitational force V4> are Gaussian variables. We can then apply the 
theorem stated, for instance, in Refs. [20, 28], which ensures that the conditional probability 
of zero-mean Gaussian variables Xa and Xb is itself a Gaussian variable with mean 

v (x B ® Xa) 

X A ) = ) - (rX A (4.6) 

7 (X A (8) X A ) 

and covariance matrix 

(x B ® X A ) 

(x a ® X A ) 


X A <8> Xb)- 


(4.7) 


C(X B ,X B ) = (. X B ®X B ) 



In our case, we identify Xa with X5r and X B with V4>/j. As explained above, the gravita¬ 
tional force acting on the halo centers is smoothed to reflect the finite extent of the halos. 
Since (V4> (g> 5) = 0 because of rotational invariance, we can simply compute (V4>) p k as 
(V^rIV^r = 0).The mean shift of V4> at the peak position is given by 


V4>) = 

/ P k 


(V$R 


VdR = 


(v$R • VS R 

(ys R y 


= —a 


^o 2 (A) 

uf(A) 


V(5r 


pk 




(4.8) 


where we have highlighted the dependence of the spectral moments uj on the UV cut-off A. 
Since at the peak V<5 r is vanishing, no renormalization is needed. We will hereafter drop 
the subscript R whenever V<L and 5 appear because we will always refer to peaks, hence the 
smoothing is understood. Notice also that, while the renormalisation condition is imposed at 
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all times, the scale A is fixed at initial time and does not evolve. In reality, A, which is related 
to the smoothing scale R that define the Lagrangian extent of a halo, likely evolves with time 
(see e.g. [29]). For simplicity however, we will ignore this complication in the present work. 

We wish now to compute the average of V3>(r, t) ■ d^f(r,v,t) with itself at the peak 
position. On applying Wick’s theorem, the only non-vanishing term that remains is 


( (v<F(r,t) • d$f(r, v, t)^j ^ = d$f pk (r,v,t) ■ d#f pk (r,v,t)(v$(r,t) • V$(r,i)^ , (4.9) 

where we used the fact that ^Vj<I>(r, t) ■ V j<&(r,t)^ oc 5ij. This requires the knowledge of 

the covariance matrix for given the constraint at the peak. At the linear order the 
unconstrained correlators needed to evaluate eq. (4.7) read 


^V<f> <2> V4>^) = y CT -i( A ) !sx3, 
(v<D® V<$) = -|^o( A ) *3X3, 
(v<5® V<3?) = -|<rg(A) * 3 x 3 , 
(v$®V<$) = i<7?(A) l 3 x 3 - 


(4.10) 


The covariance matrix for V4> given the peak constraint is then equal to 


c v$,v$ = 

\ Jp k 3 


« 2 O m CT o( A ) , - 

a -n A ) - ) ^-3x3) 


aj(A) 


(4.11) 


and 

(vmt) • ^(o*)) pk = TrC (v<h, V$) = a 2 (a 2 ,(A) - . (4.12) 

We recall now the prescription (4.4): at the lowest order, [V$-3^/] contains only the operators 
O o and 0\ (see fig. 2). Hence, we learn that the renormalized force felt by the halos at the 
linear order reads 


first 


V3>(r, t) ■ dfif(r, v,t) = V<f>(r, t) + a 


07 


;( A ), 


°l( A ) 


Vd(r,t) ■ dtff(r, v, t). 


(4.13) 


Although the mean shift of the gravitational force vanishes at the peak and thus there is no 
extra force at the peak-by-peak level, the gravitational force at the peak receives a correction 
when the statistical ensemble average is taken, as exemplified by eq. (4.12). The extra effective 
force felt by the halos is purely of statistical origin and, therefore, does not violate the 
Equivalence Principle. Furthermore, the time-dependence of Zp\ is that of acrg/u 2 , that 
is, the amplitude of the resulting fc 2 -correction in Fourier space does not depend on time. 
Therefore, this gravity bias does not decay with time, in agreement with [3, 19]. Note that 
the same result eq. (4.13) was found in Refs. [19, 21, 22, 30], but our Boltzmann approach 
combined with a renormalization procedure puts it into a different perspective. 
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4.2 Renormalization at higher order with the path integral technique 

In this subsection, we will compute the renormalized operator [V3?(r, t) ■ d^f(r, v, t)] at 
the second order in perturbation theory. We will compute the correlators beyond the Gaussian 
approximation with the path integral technique (details can be found in appendix A). We 
will restrict ourselves to the four operators O o, 0 \, 0 2 , O 3 defined in fig. 2 . 

Let us give a closer look to the renormalization prescription (4.4). We begin by observing 
that the expectation values of single operator ( Oi) are not helpful, because they are all 
proportional to expectation values of vector fields and thus vanish, as in (4.5) and (4.8). 
Therefore, in order to determine the coefficients Zij of (4.2), we write down the system given 
by the expectation values of products of two operators: 

Y, ZaZjmpiOm) = (OiOj)^ . (4.14) 

l,m 

We must now specify how we define the leading and subleading orders of our computa¬ 
tion. On the left-hand side of eq. (4.14), operating at second order the correlation functions 
(OiOj) contain at most two power spectra. Therefore, at leading order we compute the cor¬ 
relators by including only their Gaussian components and obtaining one power spectrum in 
the result; the next-to-leading order includes up to two power spectra. On the right-hand 
side, the leading order of the computation is obtained by considering again only the Gaus¬ 
sian component of the fields. When computing the correlators (OiOj) p k, the linear order is 
equivalent to consider only the leading order in uaQ. 

The results for the correlators without the peak constraint are 3 


(leading order) (next-to-leading order) 


(O 0 O 0 ) = 

a 2 c 

2 a -l°ij 


(4.15) 

(OoO,) = 

« 2 e 


(4.16) 

(O 0 O 2 ) = 


a 2 

y/Ct 6ij 

(4.17) 

(O 0 O 3 ) = 


a 3 

- 1C iSii 

6 J 

(4.18) 

(OfOf) = 

1 2 r 

2 a l °ij 


(4.19) 

( 0102 ) = 


a 17 4 

— "3 ~y~ a 0^ij 

(4.20) 

(Oi0 3 ) = 


-y/C 2 <% 

(4.21) 


(4.22) 


where we are dropping the factors of d$f to simplify the notation. The correlators calculated 


3 see Appendix B for details about the calculation of these correlators, and for the definition of the finite 
quantities Ah, Ah) 
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with the peak constraint read 


(OoOo)pk 

<0O02>pk 

<0 O 03>pk 




+(rao) 


y(Wo)<% 

CK 3 


^1 




On 


0 "T 




0 "T 


(4.23) 

(4.24) 

(4.25) 


We also notice that our renormalization prescription (4.4) enforces (OiC^pk = 0 Vi. This 
corresponds to a subsystem of (4.14) (when i = 1) whose right-hand side is 0, the solution of 
which is Zu = 0 or [Vh • d$f] = 0. This is a consequence of our renormalization prescription: 
once we impose that at the position r there is a peak defined by Vh(r) = 0, the operator Vh 
after renormalization cannot be anything but zero. 

We are now free to impose the condition [Oi] = O t + .... i.e. the diagonal entries 
satisfy the condition Za = 1 for i ^ 1. Let us first discuss the system for the renormaliza¬ 
tion of the three operators {Oq, 0\, 02}. We begin by inspecting the equation (OoOo)pk = 
ZoiZoj(OiOj), which contains all (and only) the coefficients Zo%. We know that 

2 

Zoo = 1; Zoi = a—| + ©(wo), Z 02 = O{vao) ■ 


We should impose the equality (0oC\))pk = ZoiZoj(OiOj) separately between the leading 
order terms and the next-to-leading order ones. We find that 


Zoi = ol— | + O(uao), Z02 

a i 


(^oo) 


1 | 31cio 

2^ + 6 (—17<t§ + 7/Ci<t?) 


+ O ((z-'oo) 2 ) . 


We now look at the equations (Oo^ 2 )pk = ZoiZ 2 j(OiOj) and (0202) P k = Z 2 iZ 2 j(O l Oj), 
which include the unknowns Z 20 and Z 21 . These two equations, once we select the terms con¬ 
sistently with their order (for example, the connected contribution of ( O 2 O 2 ) = (hV4>5V<h) 
contains 3 power spectra and should not be put together with second order quantities) we get 


20 w-<rv;) 



Z 21 — OL 


On 


>^1 - Y a o a -i 


K - a -l a l) a l \ (.*0 ~ 0-1*1) 


(£ 1*1 - T a o y 


+ °0 


(¥) 2 


(T 6 _ fj 2 

°0 a -l 


or 


O, 


+ a -5 (ua 0 ) 


a 


1 


—<5pk-^01 


This concludes the determination of the renormalization coefficients at the next-to-leading 
order for the set of operators {O 0 , 0 1 , O 2 }, keeping up to two power spectra in the (connected) 
correlators. We could now include the operator O 3 in the discussion, since it contains the same 
number of fields <f> and derivatives acting on them as 02 - But in this way we should compute 
correlation functions as ( 0303 ) p k, whose connected part contains three power spectra, and 
thus goes beyond the perturbative order we are considering. The difference with respect to 
the operator 0 2 is that the correlator ( 0202 ) p k> with the constraint 6 = vero, gives (i'O'o) 2 
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times a correlator with one power spectrum. Now, <7 q is of the same order of a power spectrum 
in a perturbative expansion (see eq. 3.2), thus (C^C^pk is of a lower order with respect to 
( 0303 ) p k. Therefore, at the perturbative order we are considering, [V4> • d #/] does not mix 
with ((V$-V)V$) • dtf. 

We have therefore found that the next-to-leading order renormalization correction reads 


V$(r,f) • dtff(r, v, t) 


NLO 


+ 


Z'Cr 0 (A) 


+ 


L2a 0 2(A) 

3 l^o (A) 


6 [—17(7q(A) + 7/Ci(A)<Tp(A)] 



d$f(r,v,t). 

(4.26) 


Let us summarize the main results of this section, as they are the central point of this 
work: we have renormalized the composite operator [VT • d$f] imposing the peak constraint. 
Already at first order, we have found that while the mean of Vd* is zero as in the unconstrained 
case, the variance has an additional term. The composite operator therefore renormalized in 
such a way as to generate the correct variance and account for this statistical effect, in (4.13). 
We then compute the next-to-leading order effect and find the corresponding renormalization, 
(4.26). 

4.3 Non-renormalization of the DM Boltzmann equation 

It is crucial to notice that for the smooth component of DM (that is without peak 
constraint) no renormalization is needed and the corresponding Boltzmann equation is not 
altered. This is consistent with the fact that we have used standard perturbation theory 
results for DM to evaluate the statistical correlators at the peak locations. Consistency 
can be explicitly checked upon noticing that, for DM, probabilities are not conditional. For 
instance, we get at linear order 

(V$>(r,t) = f(r,v,t) =0, 

(v$(r,t) f(r,ih,t) • V$(f,t)/(r,u 2 ,t)) = (f{r,vi,t)f(r,v 2 ,t)^(v<Z>(r,t)-V${r,t)^ 

2 

= -^-<7-1 ^/(c vi,t) f(r,V2,t)^ (4.27) 

and one simply finds 

V4>(r, t) 

This remains true also at higher orders. 

5 The halo Boltzmann equation: the right-hand side 

We will now deal with the right-hand side of the Boltzmann equation, see Eqs. (2.11) and 
(2.12). It contains a force term that depends on the irreducible two-particle correlation (which, 
in turns, depends on the irreducible three-particle correlation etc.). It is the correlation en¬ 
coded by f 2c that causes the non-conservation of the one-particle phase-space distribution 


= V4>(r, t). 


(4.28) 
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function in configuration space. The physical interpretation of this effect is that the gravita¬ 
tional interactions between particles (which induce long-range correlations) lead to clustering 
and, therefore, non-conservation of the one-particle phase-space distribution. 

To calculate this extra force, we follow Ref. [28]. The procedure consists in relating phase- 
space densities to probability distributions for the mass and velocity field, which are naturally 
specified by models of cosmological structure formation. As explained in sec. 2, we can express 
the phase space density f(r,v,t) in terms of the probability distribution P(p,v,t). This 
approach can be extended to the irreducible two-particle correlation. Inspecting eqs. (2.11) 
and (2.12), we observe that f 2 C (r,v,r' ,v' ,t) is integrated over d 3 t/ and we needs only the 
velocity at one point (the superscript ' indicates that quantities are evaluated at r') 


f2c(r,v,r',t ) = j d 3 v' f 2 c(r,v,f',v',t) 


roo roo 

/ d pp d p p [P(p,v,p\t) - P{p,v,t)P(p,t)\ 

Jo Jo 

roo roc 

J dpp J dp p [P(p,p'\v,t) - P(p\v,t)P{p',t)] P(v,t) 

[{pp'\v) - ip\v){p)] P(v',t ) , 


(5.1) 


where, in the last equality, the averages are equal-time correlators. Even though one derives 
these equations assuming a single-valued velocity held in each realization, they are also valid 
when there is a distribution of velocities at each point. We simply interpret p as the total 
mass density while v is a single bulk velocity. Therefore, these equations are fully general 
[ 281 - 

Let us specialize to the case of halos. The phase-space distribution (2.7) can be written 
using Bayes’ theorem as 

roo 

f P k(r,v,t) = dppP(v\p,t)P(p,t) = p pk P(v,t\pk) = p(l + 5 pk ) J~ 3 P(tf, t|pk) (5.2) 
where, in the last passage, we have introduced the quantity 

mt) = — VS(r,t) = - / —- ^^6(r’,t), (5.3) 

a J An \r — r'| a 


and is the Jacobian describing the passage from the variable v to the variable '0. 

In a similar fashion, we can compute the irreducible two-point correlation function 

roc roc 

f^(r,v,r',t) = / d pp / dft ft [P{v,p'\p,t)P(p,t) - P(v\p,t)P(p,t)P(pf,t)\ 

Jo Jo 

roo 

= Ppk / dp’ p’ [P(v, p\ i|pk) - P(v, t\pk)P(p', f)] 

Jo 

roo 

= p pk P(v,t\pk) / dp p [P(p’,t\v, pk) - P(p’,t)\ 


= p 2 0 + Jpk) J-.~ Pi'fp, i|pk) {S(f', t ) t),pkj - (8(r', t ) 

= P f P k{r, v, t ) (d{f ‘', t) ^(r, t), pk 


(5.4) 
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In the last equality, we have substituted (S(r', t )) = 0 as this average vanishes in the absence 
of peak constraint. The corresponding force felt by the halos is therefore 


Ppk = G N pfpk(r, v, t ) / dV U(r', t) <$(/, t), pk 


r — r 


r — r 


713' 


(5.5) 


The computation of the force (5.5) is reported in Appendix C and we provide only the final 
result up to second-order 


F pk (r, v,t)=- 47 t G n p f p k(r, v, t) 




- 2 V<J (^0 




S(r,t) a\ 
a o - a l a -i a o 


X ( (6(r,t)i/)(r,t) ■ i/>(r,t)} - J d£ ^ (5^, t)6(r, t) 


2 

+ [ifi- -=^V5(r,t) 


cr, 


o 


t(r,t) 
a o “ a l a -i 




(5.6) 


— ( 5(r,t)ijj(f,t ) • VS(r, t) 


where l = {r' — r). Using the results for the cumulants reported in (B.9), (B.10), (C.6), (C.7) 
we find (notice that the last two lines in (5.6) cancel each other) 


Fpk(r, v, t) = -4vr G N p / pk (r, v, t ) 


ip{r,t) + ( ij}(r,t) - V<5(r,t) 


5(f,t) 


07 


4 _ 2 2 

°0 a l cr -l 


CTn 


10 
— ( 
21 


”r — 22 
ba-i — 7rr (J o cr -i 


• (5-7) 


Two comments are in order. Firstly, we can easily recover the DM case away from the peaks, 
either by computing directly the force from expression (5.4) with aid of (h(r / , t)ip(r, t)), or by 
eliminating the peak constraint taking <Tq, erf 3> 1. In both cases, we obtain 

4n (r, V, t ) = -47r Gn p/dm(u V, t) VKu i). (5.8) 


Secondly, neither the force on the peaks nor on the DM contributes to the momentum equation 
once the integration of the force over momenta is considered, 


J d\’ FV- ■ Tp k , dm (/r, v, t ) = 0, 


(5.9) 


in agreement with the Equivalence Principle. Technically, this cancellation follows from the 
fact that integrating out the velocities is equivalent to an integration over the configurations 
of the held ij). As the phase-space distribution is proportional to we obtain 


d 3 t/> t)if = (rfi) = < 


°i(A) 


pk 


= 0 (for peaks) 
(for DM), 


(5.10) 


where, for peaks, we made use of eq. (4.8). This equation seems to suggest that one would 
break the equivalence principle if the constraint were different from V<5 = 0, as may be the 
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case for other tracers of the large scale structure. However, isotropy and homogeneity would 
enforce a constraint on the modulus of Vd only, e.g. |Vd| < c for some constant c > 0. 
Therefore, the integral over ip would also vanish in those cases, and the Equivalence Principle 
would still be satisfied, as it should be. 


6 Discussion and conclusion 

The fine understanding of the large scale structure of our Universe is a central theme in 
cosmology. Our knowledge in the field is nowadays (and will be even more in the future) 
dominated by data. It is, therefore, essential to understand in detail the statistical properties 
of the observed galaxies. Achieving this goal automatically calls for a detailed investigation 
of the statistical properties of the DM halos, where galaxies are expected to reside. 

In this paper we have derived an effective Boltzmann equation that describes the dy¬ 
namics of the halo mean-field phase space. This equation is necessarily different from the 
Boltzmann equation of DM as halos are not statistically unbiased relatively to the coarse¬ 
grained DM density field. Through a renormalization procedure of composite operators, we 
have obtained up to second-order in perturbation theory the statistically biased gravitational 
force felt by halos. We conclude that, at this order, the Boltzmann Equation for halos reads 


dfpk 

dt 


+ v- 


dfpk 

dr 


V$ R (r,t) + a 


°l(R) 


\75 R (r,t)+ 


+ v<j 0 (R) 


1 _ 31oq(-R) _ 

2<7q(.R) 6 [—17 a§(R) + 7/Ci( J R)af(i?)] 


$R(r,t)V$ R (r,t) 


dfpk 

dv 


= 0 , ( 6 . 1 ) 


where the extra contributions to X7& R (r,t) are purely of statistical origin: they do not arise 
because halos feel a different gravitational force than dark matter, as this would violate the 
equivalence principle, but rather as a result of imposing the peak constraint on their statistics 4 . 
This should be contrasted to the DM Boltzmann equation 


df 

dt 



V $ (r, t ) 



( 6 . 2 ) 


We emphasize again that / p k is the one-particle phase space density associated with the mean- 
field distribution of halos that corresponds to a particular /, that is, a particular realization 
of the large scale structure. The biased gravitational force experienced by the halo centers 
imprints a signature, among others, in the peculiar velocities of virialized halos. This generates 
the first-order statistical velocity bias which has been discussed in [3, 19, 21] and measured 
in [19, 31]. However, tidal forces etc. will also be affected. Our approach can in principle be 
applied to work out these corrections at any order in perturbation theory. 


4 An alternative, perhaps more intuitive, interpretation of this result is the following: instead of modifying 
the Boltzmann equation for halos with an effective force field, one may interpret the results of equations (4.13) 
and (4.26) as the corrections to the probability distribution function (PDF) of the force acting on halos. In 
this picture, the Boltzmann equation for halos (6.1) would look exactly as the one of DM, equation (6.2), but 
with a conditional PDF of the form of the terms in bracket in equation (6.1). 
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To give intuition about the magnitude of the statistical correction to the gravitational 
force in eq. (6.1), we plot in fig. 3 the amplitude of the second and third terms in parenthesis 


UR) 

UR) 


a 


goCg) 

UR)' 


isa 0 (R) 


+ 


31 UR) 


2oq(R) 6 [-1 7<rg(R) + 7XIi(i2)<r?(i2)] 


(6.3) 

(6.4) 


as a function of the Lagrangian scale R = 1 /A assuming the biasing is described by the peak 
constraint. The magnitude of the second-order statistical correction also rises with R, in 



R (Mpc) 

Figure 3. (R) and a&iR) plotted here as a function of the Lagrangian scale R = 1/A. Here we 

normalize £ 2 ( R ) with a factor of a which relates S with V/> in (6.1). 

agreement with the expectation that these effects should increase with halo mass. 

Measurements from N-body simulations are challenging since the effect must vanish in 
the limit k —>• 0. Furthermore, the sparsity of massive halos could bias the measurements 
when the purpose is to define volume-weighted statistics (see e.g. [32] for recent discussions). 
Note, however, that number-weighted statistics are most suited to extract these effects as 
they are the quantities being observed. 

There is some confusion in the literature regarding the nature of this effect. A physical 
bias will immediately arise as soon as there are at least two different fluids since, in this case, 
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velocities do not need to share the same value. This bias is effective on an object-by-object 
basis. This is obviously the case of e.g. baryons and DM. Furthermore, this will also be 
the case of the halos and DM since, owing to the finite halo size, the force acting on the 
halo center-of-mass, V<L^, is different from the force acting on DM particles, VT. Such a 
coarse-graining procedure can introduce additional terms in the fluid equations, as was already 
pointed out in [33]. This possibility was recently reconsidered in [34], If we ignore in Eq.(6.1) 
all the statistical corrections but the term then our results are identical to his, except 

for the time constancy of our R. As shown in [34], a time-dependent filter W(f/R,t ) would 
yield extra terms in the moments of the Boltzmann equation owing to the time derivative 
df/dt. 

However, [34] did not include the statistical corrections, which are the focus of this work. 
In fact, our results would hold even in the extreme case where no extra filtering is applied to 
the macrostate. This illustrates the fact that the statistical effects we are discussing are not 
caused by the finite extent of the halos. Regarding the existence of these statistical effects, 
we note that, while a time-dependent window W(f/R,t ) mimicking the collapse of halos can 
induce /^-dependent contributions at low redshift, corrections are small at high redshift [34], 
Therefore, this cannot be the explanation for the suppression of the proto-halo velocity power 
spectrum measured in [19, 31]. 

To conclude, note that halo number conservation implies that we follow the collapse of 
Lagrangian “peak-patch”, as coined by [35], into parent or isolated DM halos. Therefore, we 
do not explicitly track the merging history of these host halos. Although one could naively 
apply the extended Press-Schechter (EPS) approach [36] to determine the merger history, it 
would be very interesting to assess whether the statistical effects considered here also affect 
merger rates etc. However, it is not obvious how to do this within our current formulation, 
especially since the smoothing scale R is linked to the parent halo mass and, thus, is not free 
unlike in EPS theory. Our assumption that the density peak retain their initial critical height 
(i.e. v does not depend on time) is another limitation of our approach. A better treatment 
would consistently track the time evolution of both R and v such that the mass enclosed 
within R always equates the host halo mass. We leave all this for future work. 
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A Conditional probabilities 

In this Appendix we present the techniques we use to compute conditional probabilities 
of initially Gaussian fields which develope non-Gaussian fluctuations as time evolves. This 
method is inspired by the path integral approach and a similar application can be found in 

[ 371 - 

Let us suppose there are n stochastic variables X *, with i = l,--- ,n. We define the 
connected two-point correlators as 

Mi? = (AjXj) c (A.l) 
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and the connected three-point correlators as 


l^ijk — {XiXjXk) c 


(A.2) 


Thanks to Bayes’ theorem, we can define the conditional probability that X\ assumes a value 
given that the other variables assume defined values as 


P(X 1 \X 2 ,...,X n ) = 


P(x u x 2 ,...,x n ) 


where 


P(x 2 ,...,x n ) ’ 
P(X 2 ,..., X n ) = [ dXx P(x lj x 2 , ...,X n ). 


(A.3) 


(A.4) 


The basic quantity to compute is therefore P(Xi,X 2 , • • • , X n ). Following Ref. [37], it can be 
expressed in the following way 


P(X 1 , X 2 , ...,X n )= I V X e iX ‘ Xi 


where 


VX = 


2 vr 


dA n 

~2V 


Using the fact that 


\e iXX = -id x e tM , 


i\X 


(A.5) 


(A.6) 


(A.7) 


perturbing for small three-point correlators and defining with subscript g the Gaussian coun¬ 
terparts, we find 


P{x i|X 2 ,..., X n ) = Pg(Xi|X 2 , ...,X n ) 

6-Pg(A 2 ,..., X n 
P e (X 1} X 2 ,...,X n ) 


6P g (X ^.., X n ) didj9k '' • ’ 


(A.8) 


+ 


6P|(A 2 ,..., X n ) 

The Gaussian probability can be expressed as 


Ahjfc J dXi didjdkP g (Xi, X 2 ,..., X n ). 


P s (X 1 ,X 2 ,...,X n ) = 


1 


(2vr) Ti / 2 (Det C) 1/2 
where C is the covariance matrix. We now write 

Hmqrdmdqd r = 

T “&jl\q r d\dqd r 
+ ?>jlii r dld r 
T p>mqrdmdqdr i 


-pXiC^Xj 

g 2 b tj J 


(A.9) 


(A.10) 
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where the hat in the correlator ji mq r means that we are taking indices different from 1 . 
Using equation (A. 8) we get 


(X 1 \X 2 ,...,X n ) = 


that can be simplified 


(Xi|X 2 ,, X n ) 


M 111 


J dX 1 X 1 P(X 1 |X 2 ,...,X n ) 

(X 1 \X 2 ,...,X n ) g 

J dX 1 X 1 dfP g {X 1 ,X 2 ,...,X n ) 

J dX, X] d 1 d q d r P g (X 1 ,X 2 , ...,X n ) 

J dx, x 1 tfd r p g (x l ,x 2: , ...,x n ) 

J dx, X! d m d q d r p g (x u x 2 ,x n ) (A - n) 


+ 

+ 

+ 

+ 


6Pg(X 2 , .. . 

,x n ) 

dfllqr 


6P g (X 2 ,. . . 

,x n ) 

CO 


6P g (X 2 ,... 

,X„) 

I^mqr 


6Pg(X 2 ,... 

,*n) 

(XilXa,... 

j x n ) g 

6P g (X 2 ,... 

,x n ) 

(X x x 2 ,. .. 

j x n ) g 

6P g (X 2 , ... 

,x n ) 

(X x x 2 , ... 

j x n ) g 

6Pg(X 2 , ... 

,x n ) 

(Xr X 2 , ... 

) A^n)g 

6P g (X 2 , ... 

,x n ) 


Mm J dXr df P g (X 1 ,X 2 ,..., X n ) 

3Al qr J dXi didqdr Pg(X\. X 2 , . . . , X n ) 
3/illr J dXi dfd r Pg(X l5 X 2 ,..., X n ) 

firriqr / dXi An Oq () r Pg(Xi, X 2 , . . . , X n j. 


m 


= 1 dX 1 X 1 P(X 1 |X 2 ,...,X tl ) 

= (X, |X 2 , . . . , Xni) g 

+ 0 

+ 2P g (Xv".,X ri /' ?5rPg(X2 ’ • • • 5 Xn) 

+ 0 

- 6 p (x™Z.,X n ) dmdqdr S dXx x > p g( x >’ x 2, 
+ 0 
+ 0 
+ 0 

+ ‘ ’ X : )g Vm q rdmd q d r P g (X 2 , X n ). 


,Xn) 


(A.12) 


6P g (X 2 ,..., X n ) 
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Similarly, 


{Xf\x 2 , .. . , X n ) = j dX, XlP{X l \X 2 ,...,X n ) 
= (Xf\X 2 ,...,X n ) g 


+ ——^—— dqdr / dXl x 1 P g (X u X 2 , ..., X n ) 

ig(A2 , . . • , sin) J 

- , - llr - rd r PJX 2 ,..., X n ) 

P g (X 2 ,...,X n ) sy 2 ’ ’ n> 

^ mqr dmdgdr J dX i Xl P g (X,,X 2 ,..., X n ) 


+ 


6 P g (X 2 ,...,X n ) 
(Xl\X 2 ,...,X n ) g 
6 P g (X 2 ,...,X n ) 


fi'mqrdm d q d r P g (X 2 ,...,X n ). 


(A.13) 


B Cumulants at second-order 

In Appendix A we have showed how to calculate the average of a stochastic variable 
which is conditioned by a number of other variables, taking into account its non-Gaussian 
nature. 

In this Appendix we apply such a technique to compute the average of Vd> on the peak 
up to second-order. This amounts to using (A. 12), which requires the knowledge of the the 
cumulants [i mqr = {X m X q X r ), where the Xj range among the seven Gaussian variables <5, 
Vjd), and Vj<5 (i = 1, 2, 3). Here we are omitting the dependence on r and t as the variables 
are all evaluated at the same point. Because of rotational invariance, the only non vanishing 
cumulants are 


(<5(r, t) 5{r , t) S(f, t)), (5(r , t ) V*$(r, t) Vj$(r, t)), 

(5(r, (5(r, t)Vi6(r,t)VjS(r,t)). 


(B.l) 


We show the explicit computation of one of them, the others are calculated in a similar 
fashion. Let us take for example (<5Vj<I> Xj5). Going to Fourier space we have 




d 3 h d 3 k 2 d 3 k 3 d 3 g 

( 27 t ) 3 ( 27 t ) 3 ( 27 t ) 3 ( 27 t ) 3 


(5(q) 5(h - q) d(k 2 ) 5(k 3 ))F 2 (q, h - q)+ 


+ eye. 


ia ) {-ik 3j ) 
k 2 


= a 


— a 


— a 


d 3 L ’2 d 3 ^ 

( 27 r ) 3 ( 27 t ) 3 

f d 3 Aq d 3 /cs 
/ ( 27 r ) 3 ( 27 t ) 3 

f d 3 h d 3 k 2 
I ( 27 r ) 3 ( 27 t ) 3 


P(k 2 )P(k 3 )2F 2 (k 2 ,k 3 ) 
P(h) P(k 3 )2F 2 (k 1 ,k 3 ) 
P(h) P(k 2 )2F 2 (k 1 ,k 2 ) 


h jk 3j 

k 2 

( (hj + k 3i )k 3j 

V \h + k 3 \z 

hj (hj +k 2j ) 

k\ 


where we have used the usual standard perturbation theory kernel 


5 , 1 9 i ■ 92 (Qi 


^2(91,92) = = + x 

7 2 qiq- 2 


92 


-1-I + « 

92 91 


2 / f/i • 92 


9192 


+ 


(B.2) 


(B.3) 
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and the Poisson equation 


V 2 d> = a6. 


Let us define A by 
By taking the trace we get 


(SViQVjS) = A 5ij. 


A 


a 

3 


^-^klklnk^nki) J dn2F 2 (kiM) 


1 (fci + k 2 ) • k 2 
\h + k 2 \ 2 


(B.4) 

(B.5) 


(B.6) 


where /i = cos 9, being 9 the angle between k\ and k 2 . We can now integrate over n and 


obtain 

A = -a^ J 2 * 22 ^ 2 *1 fc 2 p ( k i) p {k 2 ) = -a (B.7) 

A similar calculation of the other cumulants gives 

Q4 

(SSS) = y a 0 4 , (B. 8 ) 

a 2 

= —K-iSij, (B.9) 

(SV&VjS) = (B.10) 

ry ry 

(V^V^V/V^) = — — (<5 Vjd> Vj<L) = — — /Ci dij , (B.ll) 

a 2 

(Vi<JVi$VzVj$) = -—K*8ij. (B.12) 

Here we have defined the quantities 


/Ci = 


’dfcidfc 2 

47 T 4 


P(k 1 )P(k 2 )x 


’5 

2kik 2 

k\ + kfj 

^3 kf — 14 k\k 2 + 3 k 2 

+ 3 {kl-klY (logtM) 

’ 



168fcf/c| 





/C 2 = 


’dliqd^ 


+ 


47T 4 

(^1 + fci ) 


P(ki)P(k 2 ) 


(B.13) 


2 1.2 


rfcffc 


1 ^ 2 ' 


5 

2k\k 2 

[ k i + fc ') 

^3fcf — lAkfk 2 + 3 k 2 

+ 3 (.kl-klY (log'll) 

336 A;?fco 


These quantities are finite: the integrand is finite for k± —> k 2 . and for k\ or k 2 —>• 0 the 
power spectra lead to a finite result. 

We are now ready to compute the average of V4> up to second-order. We use (A.12) 
which contains the Gaussian probabilities, encoded by the covariance matrix for the variables 
(5, V<h, V<5) (we understand the symmetric components of the covariance matrix) 
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By using (B.5) for all the cumulants we find 
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In a similar fashion, we can compute the covariance 
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(B.15) 


= (B.16) 


C The peak force 

In this section we report the intermediate results needed to compute the conditioned expected 
value (8(r',t) pk) appearing in (5.5). 

The formula needed to compute this conditioned probability is again (A. 12). In this case, 
differently from what discussed in the previous section, also 8{f') (or, shortly, S') appears 
among the variables and hence in the covariance matrix. In order to write down the new co- 
variance matrix corresponding to (B.14) we order the variables as {<5(r / ), S(r), V<h(r), V<5(r)} 
and we denote £ = (f' — r). so that the covariance matrix reads 


/ 


C = 




On 


d£(£) £\ 

d£ I 

0 


1 9 

r -‘ l3 


where we have dehned 


CM = J 

ij{r) = / 


d 3 k 

(27r) 3 

d 3 k 

(27r) 3 


1 9 

r° ls 

1 9 

r ri 3 


P(k)j 0 (kr) 
ji(kr) 


(C.l) 


p(ky- 


kr 


(C.2) 

(C.3) 


and j n (x ) are the spherical Bessel functions. The following formulae for their integrals over 
the solid angle will be needed: 
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We can now write down the final result after the application of (A. 12). We mark in blue 
the terms that vanish with the subsequent integration J d 3 * d£. We also denote the cumulants 
with a compact notation, so that for example Ts'jfys = ' V5(r)). 
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After performing the integral prescribed in (5.5), we obtain the result (5.6). 

In (5.6) there are two cumulants containing S(r'), ff s , s ^ and We report here the 

result of the integrals 
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where /C i is defined in 

(B.13). 
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